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A nonlinear response theory is developed and applied to electrostatic interactions between spher- 
ical macroions, screened by surrounding microions, in charge-stabilized colloidal suspensions. The 
theory describes leading-order nonlinear response of the microions (counterions, salt ions) to the 
electrostatic potential of the macroions and predicts microion-induced effective many-body interac- 
tions between macroions. A linear response approximation [Phys. Rev. E 62, 3855 (2000)] yields an 
effective pair potential of screened-Coulomb (Yukawa) form, as well as a one-body volume energy, 
which contributes to the free energy. Nonlinear response generates effective many-body interactions 
and essential corrections to both the effective pair potential and the volume energy. By adopting a 
random-phase approximation (RPA) for the response functions, and thus neglecting microion cor- 
relations, practical expressions are derived for the effective pair and triplet potentials and for the 
volume energy. Nonlinear screening is found to weaken repulsive pair interactions, induce attractive 
triplet interactions, and modify the volume energy. Numerical results for monovalent microions are 
in good agreement with available ab initio simulation data and demonstrate that nonlinear effects 
grow with increasing macroion charge and concentration and with decreasing salt concentration. In 
the dilute limit of zero macroion concentration, leading-order nonlinear corrections vanish. Finally, 
it is shown that nonlinear response theory, when combined with the RPA, is formally equivalent to 
the mean-field Poisson-Boltzmann theory and that the linear response approximation corresponds, 
' within integral-equation theory, to a linearized hypernetted-chain closure. 
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I. INTRODUCTION 

> 

\Q ■ Electrostatic interactions between charged macromolecules dispersed in an electrolyte solvent have at- 

tracted sustained cross-disciplinary interest because of their fundamental role in governing the physical 
properties of colloidal suspensions 0,0,0], polyelectrolyte solutions 0,0, and many biological systems. 
Colloids (nm-/im-sized particles) and polyelectrolytes (charged polymers) can acquire charge in solution 
through dissociation of counterions. Familiar examples of charged colloids are latex or silica microspheres, 
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' clay platelets, and ionic micelles suspended in water. Common polyelectrolytes are polyacrylic acid, found 

in gels and rheology modifiers, and biopolymers (e.g., DNA, proteins, starches) in aqueous solution. In all 
of these systems, bare Coulomb interactions between charged macromolecules (macroions) are screened by 
counterions and salt ions (microions). This paper formulates a general response theory of microion screen- 
ing and applies the theory to microion-induced effective pair and many-body interactions between colloidal 
' O ' macroions in suspension. 

In recent years, experimental reports of apparent attractions between like-charged macroions have focused 
attention on electrostatic interactions in strongly charged, deionized suspensions. Observations of anomalous 
thermodynamic behavior, such as bulk phase separation 0,0,0,13 an d metastable crystallites 0], an d 
direct measurements of attractive interactions between confined macroions have motivated a variety 

of theories and computer simulation methods, which have been recently reviewed |l3lll4lll5lll6|. 

A variety of simulation methods have been applied to explore effective interparticle interactions and phase 
behavior in charged colloids. Standard molecular dynamics 0] and Monte Carlo [Tol Il9l l20j algorithms 
have been used to investigate crystallization in effective one-component pairwise-interactin g sy stems, while 
powerful ab initio (classical Car-Parrinello) 0, and multi-component Monte Carlo \'2'i l24l I2M I2r3 | 
techniques have modeled effective interactions and, to a lesser degree, phase behavior. 
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Theoretical approaches can be broadly distinguished by the extent to which they include correlations 
between microions. Many approaches are founded on the Poisson-Boltzmann (PB) equation for the electro- 
static potential, which is derived from mean-held approximations that neglect microion correlations. The 
classic theory of Derjaguin and Landau [27| and Verwey and Overbeek [28| (DLVO), based on a linearization 
of the PB equation, predicts that widely separated macroions interact via a purely repulsive effective elec- 
trostatic pair potential of screened-Coulomb (Yukawa) form. Similar effective interactions have been derived 
within the frameworks of density- functional (DF) theory [12, HE EE El i response theory [IE EE and 
extended Debye-Huckel (DH) theories [33, [3a] . These more recent approaches also clarify the importance of 
a one-body volume energy [IE EE El HE HE US HE El] , which contributes a state-dependent term to the 
free energy and thus can influence thermodynamic behavior. 

Microion correlations, while often weak for monovalent microions, generally cannot be ignored in the case 
of multivalent microions, as emphasized in several recent studies of charged colloids and polyelectrolytes 
EE HE HE EE EH- Microion correlations can induce short-range attractions, which have been linked to 
condensation of DNA and other polyelectrolytes [IE HE HE EE] • Another wide class of theories that include 
some microion correlations is the class of integral-equation theories 

HE H EE EE El EE El , which predict 
multi-component correlation functions from the Ornstein-Zernike relation combined with various closures. 

Many theoretical approaches rely, in practice, on some manner of linear approximation. DLVO theory 
and linearized PB cell models [Hoi l5ll l52| are based on the linearized PB equation. The DF [29l . l30l.l3ll] and 
response theory [32ll33l.l34j approaches involve truncating expansions of free energy functionals or of microion 
density profiles. While linear approximations can be justified under a wide range of conditions, their validity 
may be questioned for concentrated suspensions of highly charged macroions at low salt concentrations (ionic 
strengths) - precisely those conditions under which anomalous phase behavior has been reported. On the 
other hand, many nonlinear theories, such as the full PB theory [EE l54l H3 and integral-equation theories, 
present severe computational challenges. In fact, the nonlinear PB equation usually yields to numerical 
solution only within cell models with simplified boundary conditions. 

The main purpose of the present paper is to extend response theory to include leading-order nonlinear mi- 
croion screening and to apply the extended theory to systematically test the linear-screening approximation. 
This extension necessarily entails three-body effective interactions between macroions and corrections at the 
pair and one-body levels, for which computationally practical expressions are derived. The predicted effec- 
tive interactions could, in future studies, be input directly into statistical mechanical theories or simulations 
to study influences of nonlinear screening on phase equilibria and other phenomena. 

The key qualitative conclusion of the paper is that nonlinear effects can significantly modify effective 
interactions, becoming increasingly important with increasing macroion charge and concentration and with 
decreasing salt concentration. Numerical calculations for bulk suspensions are performed to quantify parame- 
ter ranges wherein linearization is justified. Comparison is made with a similar extension of the DF approach, 
recently applied to wall-induced effective pair interactions [EE H3 an d to effective triplet interactions [58| • 

Outlining the remainder of the paper, Sec. [H] defines the model colloidal suspension; Sec. II I II develops 
a general response theory for the system; Sec. IIVI presents analytical results for leading-order nonlinear 
corrections to the effective microion-induced interactions; Sec. [V] presents numerical results, for selected 
parameters, and comparisons with predictions of linear response theory; Sec. IVII summarizes the paper; 
and finally the Appendix compares response theory with two related approaches, namely PB theory and 
integral-equation theory. 

II. MODEL 

The system of interest comprises colloidal macroions, counterions, and salt ions dispersed in a solvent 
(Fig.0. This multi-component mixture is modeled here as a collection of N m charged hard-sphere macroions, 
of valence — Z (surface charge — Ze) and radius a (diameter a = 2a), and N c point counterions of valence 
z in an electrolyte solvent of volume V at temperature T. Global charge neutrality constrains macroion 
and counterion numbers via the relation ZN m — zN c . For simplicity, we assume a symmetric electrolyte 
consisting of N s point salt ions of valence z and N s of valence —z (same valence as counterions) in a 
uniform solvent. The microions thus number N + — N c + N s positive and N— — N s negative, for a total 
of = N c + 2N S . The solvent is approximated, within the primitive model, as a dielectric continuum, 
characterized entirely by a dielectric constant e. 
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The macroion charge, which may be physically interpreted as an effective (renormalized) charge, is assumed 
to be fixed and distributed smoothly over the particle surface. Charge discreteness can be reasonably 
neglected if the distance between neighboring macroion surfaces much exceeds the typical distance between 
charge groups on a macroion surface, roughly a j\[Z ' . The assumption of point microions limits the model 
to systems with large size asymmetries. Furthermore, we neglect polarization effects, e.g., charge-induced 
dipole interactions |59ll60l l61 | . which are shorter-ranged than charge-charge interactions, and which vanish 
if solvent and macroions have the same dielectric constant (i.e., are index matched). 

III. THEORY 
A. Reduction to One Component 

The response theory of effective interactions is fundamentally based on a reduction of the multi-component 
mixture to an equivalent one-component system by integrating out the degrees of freedom of the mi- 
croions . In this reduction, the macroions are regarded as applying an "external" potential that per- 
turbs the (otherwise uniform) microion distribution. For a sufficiently weak perturbation (dilute or weakly 
charged macroions), the microions respond linearly. The linear response approximation has been discussed 
in refs. |32l l33t l34j|. Upon increasing the macroion charge or concentration, however, nonlinear microion 
response becomes increasingly important. This motivates the current extension of response theory from 
linear to nonlinear response. 

To simplify the derivation, we first consider salt-free suspensions and introduce salt ions only at the end. 
The model system is then described by a Hamiltonian H that decomposes naturally into three terms: 

H = H m ({R}) + H c ({r}) + H mc ({R}, {r}), (1) 

where {R} and {r} denote the coordinates of macroions and microions, respectively. The first term on the 
right side of Eq. JQ) is the bare macroion Hamiltonian, given by 

1 N m 

H m = H HS {{R}) + - tw(|Ri-B.j|), (2) 

where i?ns is the Hamiltonian for neutral hard spheres (the macroion hard cores) and v mm (r) = Z 2 e 2 /er, 
r > er, is the bare Coulomb pair interaction between macroions. In the primitive model, the solvent acts 
only to reduce the strength of Coulomb interactions by a factor 1/e. The second term of the Hamiltonian, 

H C = K C + - ]T Ml'i-rjl), (3) 

describes the counterions alone, having kinetic energy K c and interacting via a Coulomb pair potential 
v C c{t) = z 2 e 2 /er. The third term is the macroion-counterion interaction energy: 

H mc = J2Ylv mc (\R i -r j \), (4) 

i=i j=i 

where v mc (r) is the macroion-counterion electrostatic pair interaction: v mc (r) = Zze 2 /er, r > a. For 
impenetrable hard-core macroions, the form of u mc (r) inside the core is arbitrary and can be specified so as 
to minimize counterion penetration inside the cores (see Sec. IIV All . For later reference, we note that the 
macroion and counterion Hamiltonians [Eqs. J3J) and may be expressed in terms of Fourier components 
using the identity 

N 1 

v(\r i -r j \) = -'£v(k)[p(k)p(-V-N}, (5) 
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where v(k) is the Fourier transform of a pair potential v(r), p(k) is the Fourier transform of the appropri- 
ate (macroion or counterion) number density operator p(r) = Yli=i^( r " r «)j an d the Fourier transform 
convention is 

p(k)= f drp(r)e- lk r , (6) 



k 

The inverse transform is expressed as a summation, rather than an integral, in anticipation that charge 
neutrality will necessitate singling out the k = component for special treatment. 

At constant temperature and volume, the thermodynamic behavior of the system is governed by the 
canonical partition function, 

Z = «exp(-/3#)> c > m . ( 8 ) 

where /3 = 1/fesT and ( ) c and ( ) m denote classical traces over counterion and macroion coordinates, 
respectively. The two-component mixture of macroions and counterions can be formally mapped onto an 
equivalent one-component system of "pseudo-macroions" by performing a restricted trace over counterion 
coordinates, keeping the macroions fixed. Thus, without approximation, 

Z = (exp(-/? J ff off )) m , (9) 

where H c fi — H m + F c is the effective Hamiltonian of the equivalent one-component system and 

F c = -k B Tln (exp [-0(H e + H mc )]) c (10) 

can be interpreted as the free energy of a nonuniform gas of counterions in the presence of the fixed macroions. 
To simplify notation, we henceforth omit the subscript c from the trace over counterion coordinates: ( ) = 



B. Response Theory 

Although the one-component mapping is exact, the challenge now shifts to approximating the counte- 
rion free energy F c . Progress can be made by re gard ing the macroions as an "external" potential for the 
counterions and invoking perturbation theory |32|, |63| : 

F C = F + f dX (H mc ) X) (11) 
Jo 

where Fq = — fc^Tln (exp(—/3H c )) is the reference free energy of the counterions in the presence of neutral 
(hard-core) macroions (the counterions then being unperturbed, except for exclusion from the macroion 
cores), ( ) A denotes a trace over coordinates of the counterions in the presence of macroions charged to 
a fraction A of their full charge, and the A-integral adiabatically charges the macroions from neutral to 
fully charged. Although each term on the right side of Eq. is infinite, the infinities cancel to yield a 
finite counterion free energy. When the macroions are uncharged, the surrounding "sea" of counterions has 
uniform density, neglecting any confinement-induced structure, which is reasonable for typical counterion 
concentrations in colloidal suspensions (see below). As the macroion charge is "turned on," the counterions 
respond, redistributing themselves to form a double layer (surface charge plus neighboring counterions) 
surrounding each macroion. 

In practice, it proves convenient to convert Fq to the free energy of a classical one-component plasma (OCP) 
by adding and subtracting the energy of a uniform compensating negative background. The background 
energy can be expressed as = —^N c n c v cc (0), where n c is the average density of counterions in the volume 
unoccupied from the macroion cores. Note that the infinite background energy formally cancels the infinities 
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on the right side of Eq. (|llfl . Because the counterions are excluded (with the background) from the hard 
macroion cores, the OCP has average density n c = N c /V, where 77 = ^{N m /V)a 3 is the macroion volume 
fraction and V = V(l — 7/) is the free volume. Thus, 

F c = F OCP + I dA (H mc ) x - E b , (12) 







where Fqcp = i*b + E b is the free energy of the OCP in the presence of neutral, but volume-excluding, hard 
spheres - what might be loosely called a "Swiss cheese" OCP. 

In terms of Fourier components, the macroion-counterion interaction can be expressed as 

(H mc ) x = yj^2v mc (k)p m (k) (p c (-k)) x . (13) 



k 



Evidently, {H mc ) x depends through (/5 c (k)) A upon the response of the counterions to the macroion charge 
density. Note, however, that there is no response for k = 0, since /5 C (0) = j Ay p c {r) = N c , which is fixed by 
charge neutrality for a given macroion concentration. Taking this k — > limit into account, and subtracting 
the background energy, Eq. I|13f) becomes 



{H mc ) x - E b = yj ^2 Vnic(k)p m (k) (p c (-k)) A + n c lim 



N m V mc {k) + -yWcclfe) 



(14) 



To proceed further, we apply a perturbative approximation for the macroion-induced counterion density, 
adapting a standard approach from the theory of metals [63, 0, |6^, |6(j . Defining the macroion external 
potential </> ext (r) by 

ze(f> ext (r) = / dr'u mc (|r - r'|)p m (r'), (15) 



the ensemble-averaged induced counterion density may be expanded in a functional Taylor series around 
</>ext(r) = 0| in powers of the dimensionless potential u(r) = — /3ze</> cxt (r): 

(p c (r)) = p + jr ^ / dn ■•• / dr n G("+ 1 )(r-r 1 ,... ! r-r n Hr 1 )---u(r n ). (16) 

n=l U ' J 

Here po is a constant, chosen below to ensure charge neutrality, and the coefficients 

G(,l+1) (r - rx, • • , r - r„) = lim f /" {Pc f ) (17) 

u^o \du(ri) ■ ■ ■ ou{r n ) J 

are the (n + l)-particle density correlation functions [H3 of the unperturbed (uniform) OCP. The correlation 
functions are, in turn, proportional to response functions (see below). To give a physical interpretation to 
Eq. (|16|) . the counterion density induced at any point r is the net response to macroion-generated external 
potentials, applied at sets of points {ri, . . . , r„}, and propagated through the OCP via multi-particle density 
correlations. Fourier transforming Eq. I|l(jf> . we obtain (for k 7^ 0) 

(p c (k)) = G< 2 )(fc)«(k) + ^^G( 3 '(k',k-k')fi(k')«(k-k') 

+ ^72 E G (4) (k',k",k-k'-k"Hk'Hk"Hk-k'-k") + ---. (18) 



k',k" 



The coefficients & n \ which are Fourier transforms of G^ n \ are related to the n-particle static structure 
factors of the uniform OCP via & n > — n c S^ n \ where the static structure factors are explicitly defined 
by H 

S( 2 »WES(fc),i(p c (k)p c (-k)) (19) 



G 



and 

S< n >(ki, • • ■ ,k n _!) = -L (f>M ■ ■ ■ /Sc(kn-i)p c (-ki - ••• - K-i)) , n>3. (20) 

Substituting u(k) — —j3v mc (k)p m (k) [from Eq. (|15|l ] into Eq. (|18(l . the induced counterion density can be 
expressed in the equivalent form 



(p c (k)) = X (fcK lc (fc) / 5 m (k) + i 7 ^x'(k',k-k')« mc (fc')i) m c(|k-k'|) 

k' 

x p m (k')p m (k-k') + ---, fc^O, (21) 



where 



and 



X (k) = -f3n c S(k) (22) 



X '(k', k - k') = (p 2 n c /2)S^(k', k - k') (23) 

are, respectively, the linear and the first nonlinear response function of the uniform OCP. The first term 
on the right side of Eq. Ij21|l represents the linear response approximation - linear in p m (k) - while the 
higher-order terms generate, as shown below, nonlinear corrections to both the counterion density and the 
effective interactions. Finally, since the amplitude of v mc (k) is proportional to the macroion charge, then 

(/5 c (k)) A = Ax(fc)v mc (fc)/3 m (k) + ^^ x '(k',k-k')i) mc (fc')^ c (|k-k'|) 

k' 

x p m (k')p m (k-k') + ---, fc^O. (24) 

Note that the response functions describe the response of the fully-charged OCP, and so do not depend on 
the coupling constant A. 

C. Effective Interactions 

We are now positioned to derive formal expressions for the effective interactions. Substituting Eq. (|24|l 
into Eq. (|14|l . the latter into Eq. (|12(l . and integrating over A, we obtain the counterion free energy to third 
order in the macroion density: 



F c = F OCP + n c lim 



N m v mc {k) + -^-Vccik) 



-^T, ^2x(k) [v m c(k)} 2 p m (k)p m (-k) 

+ iEE *'( k '' ~ k - k')v mc (k)v mc (k')v mc (\k + k'\)p m (k)p m (k')p m (-k - k'). (25) 
W 

Evidently, the linear and first nonlinear response terms in the expansion of {p c (k)) generate terms in F c 
that are, respectively, quadratic and cubic in p m (k). These terms can be related to effective pair and triplet 
interactions between macroions. To this end, we first identify 

&k) = X(k)[v mc (k)] 2 (26) 

as the counterion-induced macroion-macroion pair interaction in the linear response approximation |82t 
l33l 01 . In passing, we note that Eq. I|2t)|) is similar in structure and physical interpretation to induced 
interactions recently derived from a coarse-grained HNC theory |49| and from a cumulant expansion of the 
counterion partition function |68j| . Combining the induced interaction with the bare Coulomb interaction 
yields the linear-response prediction for the total effective pair interaction: 

4 2 \k)=v mm (k)+vH(k). (27) 
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Now the term on the right side of Eq. I|25|) that is second-order in p m (k) may be manipulated using the 
identity [from Eq. JSJ] 

Nm 1 N 2 

E v inl(\^i - R .l) = 777 E €d(*0/Uk)/U-k) + lim fi« (fc) - 7V m ^(0). (28) 
Similarly, identifying 

«$(k,k') = 2x'(k',-k-k')i) mc (fc)i) mc (fc')«m C (|k + k / |) (29) 
as an effective three-body interaction, arising from nonlinear counterion response, and using the identity 

N m 

£ ^(Ri-R^Ri-B*) = ^^£<^k,kO[/Wk)MkO/W-k-k') 

ijt:jjik=l k k' 

- 3 / 5 m (k)/5 m (-k) + 2iV m ], (30) 
the final (third-order) term in Eq. (|25|l may be rewritten as 

3^2 E E °eff(k, k')p m (k)p m (k')/5 m (-k -k')-^E *cff (k, 0)p ro (k)p ro (-k) 

k k' k 

, N m 

= 31 E ^ ff ) (R,-R,,R l -R^ + 2^EE^ 3 ff ) (k,k')p m (k)p„ l (-k) 

' i#j"#fc k k' 

- ^££«2fr*)-^£«3(M)Mk)AnH0. ( 31 ) 

k k' k 

Combining Eqs. (|25|l and (|3 ip . and again invoking the identity in Eq. the effective Hamiltonian acquires 
the following physically intuitive structure: 

, N m N m 
H cS = H m + - £ «eff(l R i- R il)+ 37 E «cff( R 4 - R ^ R l - R fc)+^ ( 32 ) 

(2) (3) 

where ir ff (r) and u^ ff (r,r') are, respectively, the counterion-induced effective pair and triplet interactions 
in real space and E is a one-body volume energy. In Eq. I|32fl . the effective triplet interaction is the Fourier 
transform of Eq. I|29|) , while the effective pair interaction is the transform of 

v$(k) = 4 2) (k) + Av$(k), (33) 

where 

- T^E^k') ~ ^«2(M) (34) 

k' 

is the first nonlinear correction to the linear response approximation. Note that the second term on the 
right side of Eq. I|34|l can be traced back to the requirement of charge neutrality, which necessitated special 
treatment of the k = term in Eq. I|13l) . 

The volume energy E - a natural by-product of reduction to an equivalent one-component system - 
has no explicit dependence on macroion coordinates. Collecting terms that are independent of macroion 
coordinates, the volume energy takes the form 

E = E + AE, (35) 

where 



E = Fqcp + -y^(O) + N m n c lim 



v mc (k) - ^-v£d( k ) + 2^c(fc) 



(36) 
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is the linear response approximation 33, 34] and 



AE 



6V K 



5>$(k,k')-^£ 



^(k,0) 



k,k' 



is the first nonlinear correction. On the right side of Eq. (|36[) . the second term represents the 
macroion with its own counterions. The terms in square brackets on the right side of Eq. i|36|) 
term on the right side of Eq. I|37(l originate again from the requirement of charge neutrality, 
that nonlinear counterion response generates not only effective many-body interactions, but 
to both the effective pair interaction and the volume energy. In fact, as is clear from Eqs. 

(2) 

the nonlinear corrections to v^g (r) and E are intimately related to many-body interactions 
volume energy depends nontrivially on the mean macroion density, and thus can contribute 
the total free energy of the system. 



(37) 



interaction of a 
and the second 
. We emphasize 
also corrections 
. (EU) and P7)l. 
Note that the 
significantly to 



D. Physical Interpretation 

While the mathematical manipulations of response theory are simpler in Fourier space, the physical in- 
terpretation of the theory is perhaps more transparent in real space. In terms of real-space functions, the 
induced pair interaction, in the linear response approximation, can be expressed [from Eq. (|26|l ] as 

w ind( r ) = f dr i J dr 2 X(l r i - r 2 |)u roc (n)Uroc(|r2 - r|), (38) 

where x(l r i — r 2 |) is the real-space linear response function, which describes the change in counterion density 
induced at point r 2 in response to an external potential applied at ri. Referring to Fig. |3 Eq. H38[l can 
be interpreted as follows. The external potential due to one macroion (centered at the origin in Fig. |2J) 
induces at point r 2 a change in counterion density J drix(|i"i — i" 2 |)v mc (ri). This induced density, which 
depends (via x) on the pair density correlation function of the intervening medium (OCP), then interacts 
with a second macroion, at displacement r from the first, giving rise to a counterion- induced pair interaction 
energy. The linear-response contribution to the volume energy (per macroion) associated with macroion- 
counterion interactions [Eq. has a closely related form 

w£d(°)= / dr i / dr2X(|n-r2|)»mc(n)t*nc(ra), (39) 



and a similar physical interpretation, except that the induced density now interacts back with the first 
macroion, generating a one-body energy. 

Proceeding from linear to nonlinear response, the effective triplet interaction can be expressed [from 
Eq. ®] as 

"i/g (r,r') = 2 J dri J dr 2 J dr 3 x'(n - r 3, "* 2 - r 3 )w mc (r 1 )i; mc (|r2 - r|)u mc (|r 3 - r'|). (40) 

Again the interpretation is clear: the external potentials due to two macroions (top two macroions in Fig. [21, 
separated by displacement r, induce a change in counterion density at point r 3 . The induced density, which 
depends via x' on triplet density correlations in the OCP, then interacts with a third macroion, at displace- 
ment r' from the first, contributing a counterion-induced three-particle interaction energy. Considering now 
the nonlinear correction to the pair interaction, and leaving aside the term arising from charge neutrality, 
the main contribution can be written [from Eq. (|34|l ] as 



,( 2 )^ - 



Ac.'jV') = 2y dri J Av 2 J dr 3 x'(ri - r 3 ,r 2 - r 3 )v mc (ri)v mc (r 2 )v mc (\r 3 - r|). (41 i 



The interpretation is analogous to that for the triplet interaction, except that the external potentials at 
points ri and r 2 are now associated with the same macroion. Finally, the nonlinear correction to the volume 
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energy [Eq. (|37|l ]. aside from the charge neutrality term, has the form 

A ' E= ~y i / dri / dr2 / dl " 3 X '^ Tl ~ r3 ' r2 _ r 3) v mc{ri)v mc (r2)v m c(r3)- ( 42 ) 

The physical meaning of AE is similar to that of Air ff (r) , except that now the density that is induced 
nonlinearly by one macroion interacts back with the same macroion, generating a nonlinear contribution to 
the one-body energy. 

E. Random Phase Approximation 

Further progress towards practical expressions for effective interactions requires specifying the OCP re- 
sponse functions. For charged colloids, the OCP is typically weakly correlated, characterized by relatively 
small coupling parameters: T = Ab/cl c "C 1, where Ab = I3z 2 e 2 /e is the Bjerrum length and a c — (3/47rn c ) 1 / 3 
is the counterion-sphere radius. For example, for macroions of valence Z — 500, volume fraction r\ = 0.01, 
and monovalent counterions suspended in salt-free water at room temperature (Ab — 0.714 nm), we find 
r ~ 0.02. For such weakly-correlated plasmas, it is reasonable - at least as regards long-range interac- 
tions - to neglect short-range correlations. We thus adopt the random phase approximation (RPA), which 
equates the two-particle direct correlation function (DCF) to its exact asymptotic limit: c^ 2 \r) — —(3v cc (r) 
or c.( 2 \k) = —Air(3z 2 e 2 /ek 2 . In neglecting short-range correlations, the RPA is formally equivalent to the 
mean-field PB theory, as shown in the Appendix. Furthermore, we ignore the influence of the macroion hard 
cores on the OCP response functions, which is reasonable for sufficiently dilute suspensions. Within the 
RPA, the OCP (two-particle) static structure factor and linear response function take the analytical forms 

S{k) = l-n c c(2)(fc) = l + K 2 /k 2 (43) 



and 



x (k) = -0n c S(k)= 1+ ^ c /k2 , (44) 



where k — W Aim c z 2 e 2 / eksT '. As will be seen below, the parameter k plays the role of the Debye screening 
constant (inverse screening length) in the counterion density profile and in the effective interactions. In the 
absence of salt, the counterions are the only screening ions. The macroions themselves, being singled out 
as sources of the external potential for the counterions, do not contribute to the density of screening ions. 
Fourier transforming Eq. 144(1. the real-space linear response function takes the form 

x (r) = -(3n c [S(r) + n c h cc (r)} , (45) 

where 

8z 2 e 2 e~ KT 

M^)=-— — (46) 
e r 

is the counterion-counterion pair correlation function 69], which has Yukawa form, with screening constant 
k. Equation l|45|l makes clear that there are two physically distinct types of counterion response: local 
response, associated with counterion self correlations, and nonlocal response, associated with counterion 
pair correlations. 

At this point, we can specify the constant po in Eq. I|16fl . Combining Eq. I|44() with the long- wavelength 
limit of the macroion-counterion interaction, v mc (k) — > 4irZze 2 /ek 2 , as fc — > 0, we have 

lim \x(k)v mc (k)] - Z/z. (47) 



Thus, the linear response term in Eq. Ijltjfl already ensures proper normalization of /o c (r), which implies that 
Po = 0. 
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Proceeding to nonlinear response, we first note that the three-particle structure factor obeys the identity 
S (3 Hk,k') = S(k)S(k')S(\k + k'\) \l + n 2 c c( 3 )(k,k')j , (48) 

where c^ 3 '(k, k') is the Fourier transform of the three-particle direct correlation function. Within the RPA, 
however, and all higher-order DCF's vanish. Thus, from Eqs. (|22() . (|23ll . and (|48|l . the first nonlinear 
response function can be expressed in Fourier space as 



X '(k, k') = -^ x (k) X (k')x(\k + k'l) (49) 

and in real space as 

X'(ri -r 2) n -r 3 ) = [ dr x(|ri - r|)x(|r 2 - r|)x(|r 3 - r|). (50) 

2n 2 c J 

Higher-order nonlinear counterion response leads to higher-order terms in the effective Hamiltonian [Eq. (|32J) ] . 
For example, the effective four-body interaction takes the form 

^ } (k, k', k") - 6x"(k', k", k k' k")v mc (k)v mc (k')v mc (k")v mc (\k + k' + k"\), (51) 

where 

X"(k,k',k") = ^!n c ^( 4 )(k,k',k") (52) 

is the next higher-order nonlinear response function and 

S (4) (k,k',k") = S(k)S(k')S(k")S(\k + k' + k"\) 

x [S(|k + k'|) + S(|k' + k"|) + S(|k + k"|)-2] (53) 

is the four-particle structure factor in the RPA. Just as effective three-body interactions are related to 
corrections at the two- and one-body levels, so four-body interactions entail corrections at the three-, two-, 
and one-body levels, which (in Fourier space) are proportional to appropriate summations of (k,k' ,k") 
over the wave- vectors k, k', and k". In principle, these higher-order corrections could be computed to further 
check for convergence of the perturbation expansion. 

IV. RESULTS 
A. Counterion Density 

A practical expression for the ensemble-averaged counterion density is now obtained by substituting the 
linear and first nonlinear RPA response functions [Eqs. (|44|) and (|49|l ] into the expansion for (p c (k)) [Eq. |JHJ]. 
The result may be expressed in the form 

(p c (k))=p c o(k)-^A 7 ^p c0 (k')p c o(k-k'), fc^O, (54) 

c k' 

where 

p c o(k) = x(k)v mc (k)p m (k), k^O, (55) 

is the Fourier transform of the linear-response counterion density and v mc (fc) is the transform of the macroion- 
counterion interaction (specified below). Inverse transforming Eq. (|54l) yields 



< Pc (r))= Pc0 (r)- J dr'x(|r-r'|)[ Pc0 (r')] 2 , 



(56) 
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where 



o0 (r)=£>(|i--Ri|) (57) 



is the real-space linear response counterion density in the presence of macroions fixed at positions R^, 
expressed as a sum of single-macroion counterion density orbitals p\(r) - the inverse transform of Pi(k) — 
X{k)v mc (k). Equivalently, 

p c0 (r)= dr'x(\r-r'\)J2v mc (\r'-Ri\). (58) 

^ i=l 

Now substitution of Eqs. H45|) and (|46[) for the real-space RPA linear response function into Eqs. (|56[l and 
(|58|l allows the linear-response counterion density profile to be expressed as 



Pco(r) = (3n c Y^ 



i=l 

and the nonlinear profile as 



; -re|r— v' 

-v mc (\r - Ri|) + J- / dr' ^- — —v mc {\r' - R,-|) 



(59) 



(Pc(r)) = , c0 (r) + JL [p c0 (r)] 2 - ^- / dr' ^ Mr0 f . (60 ) 

The last two terms on the right side of Eq. l|6l)|l are nonlinear corrections to the linear profile and can be 
physically interpreted as arising, respectively from local and nonlocal nonlinear response of counterions to 
the macroion charge. 

The linear-response counterion profile in the presence of a distribution of macroions can be obtained by 
ensemble averaging Eq. (|55|l over macroion coordinates. In Fourier space, 

</5co(k)) m = Pl (fc)< / 5(k)) m , (61) 

where ( ) m again represents a trace over macroion coordinates and (p(k)) m is the Fourier component of the 
average density of macroions. In real space, the average density of counterions around a central macroion 
can then be expressed as 

(Pco(r)) m = pi(r) + n m / dr' g mm (r')pi(\r ~ r'|), (62) 



where g mm (r) is the macroion-macroion pair distribution function. The latter function may be obtained 
from integral-equation theory or simulation, with effective interactions as input. 



B. Macroion-Counterion Interaction 



To this point, the theory makes no assumptions about the type of macroion. Practical calculations 
require specifying the macroion structure, the macroion-counterion interaction, and the corresponding single- 
macroion counterion density orbital. Henceforth, we specialize to charged hard-sphere colloidal macroions. 
A convenient strategy, proposed in ref. [3(| and adopted in refs. [33| and |34[; specifies v mc (r) inside the hard 
core (r < a) so as to minimize counterion penetration. We thus assume 



~Zze 2 
er 

-Zze 



^ r ) = { r>a ( 63 ) 



ea 



a, r < a 



and choose the parameter a appropriately. In passing, we note that v mc (r) plays a role here analogous to 
that of an empty-core pseudopotential in the pseudopotential theory of simple metals |65l l66l ] . As shown 
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in refs. |30j and [33J, at the level of linear response, penetration of counterions inside the macroion cores is 
eliminated by choosing a, = Ka/(1 + Ka). This choice yields 



AnZze 2 



e(l + na)k 2 



K 



(k) = — — — - cos(fca) 4- — sin(fca) (64) 



and 



Pl (r) = <^ ~z~$m 1 + Ka ~r~' r>a (65) 
[ 0, r < a, 

which is precisely the DLVO expression for the density of counterions around an isolated macroion 
The above choice for the parameter a ensures that the linear term and first nonlinear term of Eq. H6()(l vanish 
completely inside the macroion core. The same parametrization also allows, however, the final nonlinear term 
in Eq. (|6U|I to be nonzero inside the core, although in practice the fractional penetration is at most a few 
percent. Independent of parametrization, Eqs. (|56|l . (|57l) . and (|60|l maintain charge neutrality by preserving 
the number of counterions, since J drp^r) = Z/z and J dr%(r) = x(k = 0) = 0. 

Another artifact of the present scheme, apparent from Eq. I|57|l . is that the counterion density profile 
around a given macroion overlaps the hard cores of neighboring macroions. More general parametrizations 
of the macroion-counterion interaction than Eq. I|63|) could conceivably eliminate counterion penetration 
within all cores. An alternative strategy would incorporate excluded volume constraints directly into the 
response functions, which then would more properly describe the Swiss cheese OCP. In such a scheme, 
x(|r — r'|) would strictly vanish when either r or r' falls inside a hard core. This condition - not obeyed by 
Eqs. I)44|l and l|49(l - would enforce exclusion of counterions from all macroion cores. Nevertheless, in the 
current scheme, the extent of core overlap is minor for macroion separations that significantly exceed the 
screening length which is usually the case in practice. 



C. Effective Interactions 



The effective interactions can be expressed in real space by evaluating the respective inverse Fourier 
transforms. From Eqs. (|29() and i|49|) . the effective triplet interaction is 

u£i(ri - r 2 ,ri - r 3 ) = fdrpiQri - r|)pi(|r 2 - r|)pi(|r 3 - r|). (66) 

n c J 

Equations f2"fi|l and l(2T| , combined with Eqs. (|4*4^l and yield the linear-response prediction for the 

effective pair interaction |33| . 

(2 ) , Z 2 e 2 ( e Ka \ 2 e 



identical to the familiar DLVO screened-Coulomb potential in the limit of widely separated macroions |27ll28j |. 
while Eq. (|34fl yields the first nonlinear correction 

(0 = / dr' P i(r') Pl (\r - r'|) [ Pl (|r - r'|) - ^] . (68) 

f2l (2) (2) 

The total effective pair potential is given by ir ff (r) = Uq (r) + Av; ff (r) . Note the distinction between the 
effective pair potential, which is the interaction between a pair of macroions in a colloidal suspension of 
arbitrary concentration, and the potential of mean force, which is the interaction between an isolated pair 

(2) 

of macroions, i.e., the low-density limit of ir fi ■ (r). 

The volume energy can be expressed - by combining Eqs. (ffi>|) . JSSJl, (|37Jl . and - as the sum of 

the linear response prediction |33| . 

Z 2 e 2 k N c k B T , . 

E = F OCP -N rn — — a -£- , 69 

2e 1 + na 2 
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and the first nonlinear correction, 

AE = ~^r- (/ dr [pi{r)]3 ~ nc J dr [pi{r)]2 ) ■ (70) 

The first and second terms on the right side of Eq. i|tj9fl account, respectively, for the counterion entropy and 
the macroion-counterion electrostatic interaction energy, while Eq. I|70[) corrects the latter term for nonlinear 
response. The final terms on the right sides of Eqs. I|t)8l) - (|7l)l) originate from the charge neutrality constraint. 

D. Effect of Added Salt 

Finally, we generalize the above results to the case of nonzero salt concentration. The average number 
density (in the free volume) of salt ion pairs, n s = N s /V' , is supposed maintained by exchange of salt 
ions with a salt reservoir through a semi-permeable membrane. The total average microion density is then 
rip = n + + n_ = n c + 2n s , where n± are the average number densities of positive/negative microions. 
Following ref. (34|> the Hamiltonian generalizes to 

H = H rn + + H m+ + H m ^, (71) 

where is the Hamiltonian of the microions (counterions plus salt ions) and H m ± are the electrostatic in- 
teraction energies between macroions and positive/negative microions. The presence of positive and negative 
microion species requires a proliferation of response functions, Xij an d Xijky hh^ = and a generalization 
of Eq. gSJl to 

X++ (r) = -(3n+ [S(r) + n+h ++ (r)} , (72) 
X+- (r) = -/?n+n_/i+_(r), (73) 

X (r) = -f3n_ [5(r) + n-h (r)] , (74) 

where hij(r), i,j = ±, are the bulk microion two-particle pair correlation functions, which depend implicitly 
on n±. Generalizing Eq. 121|) . the ensemble-averaged microion number densities are given by 



(p±(k)) - ±x ± (fc)« m+ (fc)p m (k) + ^^ X ±(k / ) k-k')W(fc / )W(|k-k'|) 

k' 

x p m (k')/5 m (k-k') + --- (fc^O), (75) 



where we have exploited symmetries: v m + — — v m -> Xh = X hi X+h = x'-\ h c t c -' *° define composite 

response functions as x+ = X++ - X+-, X- = X — - X+-, X+ = X+++ - 2 X++- +X+ — , and x'- = x' 

2% / _ H + X-++- Substituting Eq. 1|75[) into the multi-component Hamiltonian [Eq. (ZQ L the macroion- 

microion interaction contribution can be expressed as 

H m+ + H m _ = ^^ X (fc)[W(fc)] 2 /5 m (k)/5 I ^-k) + -^^ X '(k',-k-k') 

k k,k' 

x u m+ (fc)u m+ (fc')wm+(|k + k'|) j o m (k) / 5 m (k') / o m (-k-k'), (76) 

where the linear and first nonlinear response functions are now redefined as the following combinations of 
Xtj and X ' ljk - 

X = X+ + X- = X++ - 2 X+- + X — (77) 

and 

X' - X'+ - X- = XV++ - 3x'++- + ix'+— - X-— (78) 
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The net effect of adding salt is to modify the previous salt-free results as follows. First, the average counterion 
density n c in the Debye screening constant k and in the linear response function [Eq. (|44ll ] must be replaced 

2 /ek B T and xO) -> -/3n,,S(k). The first 



by the total average microion density n«. Thus, 



nonlinear response function retains its original form [Eq. (|49|l ]. but with the new definition of k. The second 
modification is in the linear-response volume energy [Eq. i|69[) ]. which becomes |34l l70j 



En — F, 



plasma 



N„ 



Z 2 e 2 



2e 1 



k B T{N + -N_) 2 
2 N+ + N„ 



(79) 



where -Fpiasma is the free energy of the unperturbed microion plasma. Finally, the effective triplet interaction 
and nonlinear corrections to the effective pair interaction and volume energy are generalized as follows: 



AE 



N m kBT (n + - n_ ) 
6 n?, 



dr [pi(r)} - n M / dr [pi(r)] 



(80) 



= -fc B T 



(n+ - n-) 



dr'pi(r')pi(|r-r'|) Pl (|r - r'|) 



(81) 



,(3), 



ri - r 2 ,ri - r 3 J 



-knT 



(n+ — ri—) 



drpi(|ri - r|)pi(|r 2 - r|)pi(|r 3 - r|). 



(82) 



Equations H80|) - (|82[l [combined with Eq. (|65|l for pi(r)] are the main new results for nonlinear effective 
interactions. These expressions imply that nonlinear effects increase in strength with increasing macroion 
charge, increasing macroion concentration, and decreasing salt concentration, and that effective triplet inter- 
actions are consistently attractive. These results also imply that in the limit of zero macroion concentration 
(n c = n + — n- — ► 0), or of high salt concentration (rt M — * oo), such that (n + — — *• 0, the leading-order 
nonlinear corrections all vanish. This dilute limit follows naturally from the fact that, for a pure symmetric 

electrolyte, response functions related by symmetry are equal (x'+++ = x' ) X+h = x'-\ )j an d so, from 

Eq. (|78|l . the nonlinear response function x' is zero. This result - a consequence of charge neutrality that is 
analogous to the vanishing of the first nonlinear (quadratic) term in the expansion of the nonlinear PB equa- 
tion - may partially explain the surprisingly broad range of validity of DLVO theory for high-ionic-strength 
suspensions. Nevertheless, even in dilute suspensions at high ionic strength, higher-order nonlinear correc- 
tions do not necessarily vanish. For this reason, predictions of the first-order nonlinear theory in the dilute 
limit may deviate from numerical solutions of the full nonlinear PB equation, which include, by construction, 
nonlinear corrections to all orders. The first-order corrections are nonetheless valuable in signaling the onset 
of nonlinearity, as shown below. 

Equations (|81|l and l|82[l may be compared with related expressions derived via the density-functional 
approach. Equation l|81|l is similar in structure to an expression for a wall-induced effective pair interaction 
derived by Goulding and Hansen [Hill ( E q- ( 13 ) of Ref. [H) if one factor of pi(r) in Eq. Q81[) is replaced 
by a wall counterion density orbital. It should be noted that these authors neglected the bulk nonlinear 
correction derived here [Eq. I|81(l ]. which is justified only in the dilute limit of an isolated pair of macroions. 
Equation (|82|l is similar to an expression for the triplet interaction derived by Lowen and Allahyarov [58j |. 
differing only by a factor of (n+ — n_)/n M and by our excluded- volume correction in the definition of n. 



E. Analytical Expressions for Effective Interactions 

Quantitative predictions of nonlinear response theory are facilitated by reducing the effective interactions 
to computationally practical analytical forms. Substituting Eq. (|65(l into Eq. I|80|) . and evaluating the 
integrals, the nonlinear correction to the volume energy can be expressed as 



AE 



N m k B T 
6 



Z 2 K 3 n u 


f 1 >! 




( e Ka \ 


ftr 


V 1 + na J 


(An) 2 


V 1 + na J 



Ei(3/ca) 



(83) 
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where Ei is the exponential integral function |7lj 

POO 

Ei(x) = / du- , x>0. (84) 



oo e~ xu 



Similarly, the nonlinear correction to the effective pair potential can be rendered analytically. The key is 
expressing the integral Ix(r) = J dr' pi(r')[px(\r — r'|)] 2 in Eq. (|81|l in the form 1\ = T~ l {pi(k)p2(k)}, where 
p2{k) = J-{[pi(r)] 2 }, with T denoting the Fourier transform operator and T~ x its inverse. Substituting 

Eq. I^Sjl into Eq. ijHJl , A^g(r) reduces, after some algebra, to 



-KG 



Av$ (r) = h (r) — + / 2 (r) — + / 3 (r) — , r > a, (85) 

where 

/ x (r) = Ci [«(r - <j) + 1 - e^ "] + C 2 [Ei (/s(r - a)) + Ei (3/so) - Ei (ko)] , (86) 
/2(r) = -C 2 Ei(3K(r + a)), (87) 

and 

/ 3 (r)=C 2 [E 1 (2/c(r + o))-Ei(2K(r-a))] 1 (88) 

with 

Ci = l(n + -n_)gV/ 

6 e \ 1 + «a J 

and 

^l^W^V 
87r ze \1 + Ka J 

It is interesting to examine the asymptotic (r — * oo) behavior of the leading-order nonlinear response 
approximation to the effective pair interaction. From Eqs. I|85 |l -(|88 |l . using the inequality Ei(x) < e~ x /x, 
we find the asymptotic limit 

lim v%)(r) =G\Ke- Kr , (91) 

which exhibits a more gradual decay than the screened-Coulomb DLVO potential [Eq. I|67() ]. We empha- 
size that this result does not contradict measurements of DLVO-like interactions between isolated pairs of 
macroions in dilute suspensions fUj, since in the dilute limit C\ — > and the asymptotic behavior reduces 
to that of linear response. The physical significance of Eq. (|91|l may be limited, however, by the neglect 
of higher-order nonlinear terms and shielding by intervening macroions |72( at distances beyond the mean 
nearest-neighbor separation, r > [3/(47rn m )] 1 / 3 . 



V. NUMERICAL INVESTIGATIONS AND DISCUSSION 



To quantitatively illustrate the influence of nonlinear screening, we compute counterion density profiles and 
effective pair and triplet interactions for selected system parameters. All results presented are for the case of 
monovalent counterions (z = 1) and aqueous suspensions at room temperature (As = 0.714 nm). Figure 
compares the counterion density profile around a single macroion of diameter a = 100 nm and valences 
Z = 100 and 500 in a dilute (r) = 0.01) salt-free suspension, as predicted by linear response (DLVO) theory 
[Eq. (|58f) ] and by first-order nonlinear response theory [Eq. l|tjU|) ]. The linear-response counterion density 
Pco(r) is approximated here by a single orbital pi(r) [Eq. 1(55)1] . Evidently, nonlinear response sharpens the 
distribution of counterions around a macroion. 



16 



Figure 01 compares the linear and nonlinear response predictions for the effective pair interaction, now for 
Z = 400 (a) and 700 (b) and with small concentrations of added salt. As a check on the calculations, the 
interactions were computed both by Monte Carlo integration of Eq. (j81|) and from the analytical expressions 
[Eqs. l|85 M 9l) |) ]. the results being identical to within numerical error. Enhanced screening, resulting from 
the sharper nonlinear counterion profiles around macroions, has the general effect of weakening the pair 
interactions. For a given macroion diameter, nonlinear corrections increase in magnitude with increasing 
macroion valence and concentration and with decreasing salt concentration. Qualitatively, our predictions 
are consistent with the recent observations of Durand and Franck 73] of surprisingly short-ranged pair 
correlations in highly deionized colloidal suspensions. A quantitative comparison, however, would require 
computing the radial distribution function g(r) from our v cS (r), by means of integral-equation theory or 

( 2) 

simulation, or computing (r) from the experimental g(r) data. Note that the effective pair potential 
discussed here is distinct from the potential of mean force v m f(r), which was obtained in ref. |7^| from the 
experimentally measured g(r) via the definition v m {(r) = — fe^Tm g{r). 

In Fig. the parameters (macroion diameter, a = 652 nm, and volume fraction, i] = 0.0352) are 
chosen to compare with the experiments of Larsen and Grier |10| in which unusually long-lived metastable 
fee crystallites were observed. The macroion valence here is set to the maximum consistent with charge 
renormalization [j4|, Z* ~ O(10)(a/ Xp) ~ 5000, assuming the charge of a macroion to be reduced in bulk 
compared with its value in isolation |75j| . While the pair interaction remains repulsive, it is significantly 
weaker than the DLVO prediction over a range comparable to the macroion diameter. A weaker pair 
interaction could promote the influence of three-body attractions, as well as interactions ignored by mean- 
field theory, such as short-range counterion fluctuation- induced attractions |25l |40| . In this way, nonlinear 
response may contribute to explaining experimental evidence for apparent attractions between like-charged 
macroions. The phase behavior of highly-charged colloidal crystals will be the subject of a future study. 

The insets to Fig. 0] illustrate the extent to which the effective pair interaction may be fit by a screened- 
Coulomb (DLVO) potential. For sufficiently weak nonlinearity (Fig. v^(r) may be reasonably fit by a 
DLVO potential with the same screening constant but a lower (renormalizcd) effective charge. The tendency 
of nonlinear screening to preserve the DLVO form of potential is consistent with conclusions from PB cell 
model calculations |7J| , ab initio simulations 0, , and optical tweezer experiments [ljj , all of which 
support the bulk DLVO potential in the weakly nonlinear regime. In the more strongly nonlinear regime 
(Figs. EJ) and^]:), however, our calculations indicate that effective pair interactions may deviate significantly 

(2) 

from DLVO form. In this regime, if ir fi (r) can be fit at all by a DLVO potential, then it is only over 
a limited range and only by allowing renormalization of both the charge and screening constant. Similar 
departures from DLVO behavior with increa sing macroion-counterion coupling strength have been predicted 
by integral-equation theories HE EE HE E3 EB 113 ■ 

Although in the physically relevant range of macroion valence (Z < Z*) the predicted pair interaction is 
always purely repulsive, at higher (unphysical) valences (Z > Z*) the leading-order nonlinear theory predicts 

(2) 

that v„J (r) can develop an attractive well at sufficiently high macroion concentrations. Mathematical 
proofs |76l l77l l78j have recently shown, however, that PB theory cannot yield pair attraction - at least 
between a pair of isolated macroions. Since the RPA used here is formally equivalent to mean-field PB 
theory (see Appendix), the emergence of an attractive pair potential is best interpreted as a sign that 
higher-order nonlinear terms must then be included. 

As a quantitative test of the nonlinear response theory, we compare predictions with available data from the 
ab initio simulations of Tehver et al. jjjj . By assuming a counterion density orbital and ignoring counterion 
density fluctuations, the ab initio approach provides the most direct test of the theory. Figure|5ji presents the 
comparison for the total potential energy of interaction between a pair of macroions, of diameter a = 106 nm 
and valence Z = 200, in a cubic box of length 530 nm with periodic boundary conditions (taking into account 
image interactions) in the absence of salt. The theory is in essentially perfect agreement with simulation, 
although nonlinear effects, for these parameters, are relatively weak. Figure [SJa shows results for a higher 
valence (Z = 700), for which case simulation data are not yet available, but where nonlinear effects are more 
prominent. Further simulations of more highly charged macroions would more severely test the theory - in 
particular, convergence of the perturbation expansion - in the strongly nonlinear regime. 

To quantify the range of validity of the linear response approximation and to measure the impact of non- 
linear screening on thermodynamics, we calculate the magnitude of the leading-order nonlinear correction to 
the pair interaction Av^J(r) at the mean nearest-neighbor separation r nn , where pair interactions contribute 
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most to the potential energy. Figure HO maps out, in the space of macroion volume fraction 77 and salt concen- 

(2) 

tration c s (measured in /xM), the boundary of the region within which \Av^J(r nn )\ exceeds typical thermal 
energies, for the fee crystal structure: r nn /a = 2~ 1 / 2 (27r/3?7) 1 / 3 . For points above the boundary curves, 
\Av^(r nn )\ > 1 k B T (Fig. Hi) or 0.1 k B T (Fig. 03 ). Points on the boundary curves in Fig. correspond 
to thermodynamic states for which a stable fee crystal phase is predicted by simulations of model Yukawa 
systems [l7j ] , although these simulations do not include influences of the volume energy. With increasing Z 
and decreasing c s , the threshold r\ decreases. Thus, nonlinear screening is anticipated to increasingly influ- 
ence thermodynamics with increasing macroion charge and concentration and with decreasing ionic strength 
- just the conditions under which anomalous phase behavior has been observed 0, 0, El Ej ■ 

Moving beyond pair interactions, Fig. [7| shows the effective three-body interaction between a triplet of 
macroions arranged in an equilateral triangle for a = 100 nm and two different valences, Z — 500 and 
700. The interactions were computed numerically by Monte Carlo integration of Eq. I|82[l . The strength 
of the interaction is seen to grow rapidly with increasing macroion valence and with decreasing separation 
between macroion cores. In a concentrated suspension of highly-charged macroions, effective many-body 
interactions may become significant. In particular, as noted above, triplet attractions may well contribute 
to the surprising metastability of colloidal crystallites observed in deionized suspensions |l(J . 

To again test the theory against simulation, we compute the force on a macroion in an equilateral-triangle 
configuration of three macroions as a function of the triangle edge length. On expanding a triangle from 
edge length R — AR/2 to R + AR/2, the energy changes by 



AU = 3 



(2) 



R 



2 J 



2 J 



AR \ vWfp AR 



(92) 



Since, as the triangle expands, each of the three macroions moves a distance Ai?/y3 parallel to the total 
effective force F acting on it, the change in energy also may be expressed as AU = —3FAR/^/3. Equating 
the two expressions for AU, the total force can be written as 



F{R) = FW(R) + F<- 3 \R), 



(93) 



where 



F^(R) 



-V3 
AR 



,(2) 



R + 



AR 



,(2) 



R 



AR 



(94) 



and 



F^(R) = 



V3AR 



AR\ 



v cs \ R 2" 



(95) 



are the effective pair and triplet forces, respectively. 

To compare directly with available simulation data |2l| . we consider macroions of diameter a = 106 nm 
in a cubic box of length 1000 nm with periodic boundary conditions in the absence of salt. Over a range of 
macroion valences, we compute the sum of linear (DLVO) effective pair forces, the sum of nonlinear effective 
pair forces, the effective triplet force, and the total effective force (sum of pair and triplet forces), from 
Eqs. (|93(l - (|95|l . For a valence of Z = 200 - the only case for which simulation data were reported plj - 
the predicted total force is essentially identical to the sum of pair forces, consistent with ref. [2l|, in which 
an absence of many-body effects was concluded. For higher valences, however, three-body forces are not 
negligible. Figure |S] presents predictions of the theory for Z = 700 and 1000, demonstrating significant 
deviations of the total force from the sum of pair forces. For the case Z = 1000, which somewhat exceeds 
the charge-renormalization limit [74|. the predicted total force actually becomes attractive beyond r ~ 2a, 
although this is likely an artifact of truncating the perturbation series and thereby neglecting higher-order 
nonlinear terms. Again, further simulations could help to resolve the issue. 

Other workers have investigated many-body interactions in charged colloids. Schmitz [79^ has developed a 
theory that describes sharing of counterions between macroions, analogous to molecular chemical bonding, 
and used the theory to study the influence of many-body effects on counterion distributions and the structure 
of colloidal crystals. Sear |80|. exploring a phenomenological model, showed that interactions among particles 
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with internal degrees of freedom are, in general, non-pairwise-additive, and that tri plet interactions may be 
attractive at the same time that pair interactions are repulsive. Recently, Russ et al. |8l| solved the nonlinear 
PB equation for triplets of macroions immersed in an electrolyte. Their conclusion that three-body effects 
are always cohesive agrees qualitatively with our results, and those of ref. [H^. We note, however, that 
three-body contributions to the grand potential, calculated in ref. |8l|. are not directly comparable to three- 
body interactions in the effective Hamiltonian, calculated here and in the simulations of Tehver et al. |2l( . 
In another study, Wu et al. extracted three-body forces from Monte Carlo simulations of macroion 
triplets in equilateral configurations. These authors found attractive electrostatic three-body forces, but also 
detected a significant repulsive contribution attributable to hard-sphere collisions between macroions and 
microions, which were modelled in the simulations as charged hard spheres. Future extension of the response 
theory from point microions to hard-core microions would allow a more direct comparison with these Monte 
Carlo data. 

Influences on thermodynamic phase behavior of nonlinear corrections to both effective interactions and 
the volume energy are now being explored. It may be anticipated that these corrections will be especially 
significant for deionized suspensions of highly charged macroions. Preliminary calculations of free energies 
and pha se diagrams [83j indicate that the spinodal-instability mechanism proposed to describe phase separa- 
tion |30L |3l| remains qualitatively valid - at least under the assumption of fixed macroion charge - but that 
nonlincarity can significantly shift the phase boundaries and, in some cases, enhance the tendency toward 
phase separation. 

The role of nonlinear response and of effective many-body interactions in dense electron-ion (metallic) 
systems has long been recognized and discussed [841 Ha, |8(| H3, ■ In this context, it has been argued that 
nonlinear corrections to pair potentials and structure factors either are weak |88j or can be incorporated into 
the linear response scheme |87|. but that nonlinear corrections to the volume energy and thermodynamic 
properties (e.g., bulk modulus) are more significant |86[ . Whether the same argument applies also to colloidal 
systems is being investigated in ongoing studies of phase behavior |8fl| . 

VI. SUMMARY AND CONCLUSIONS 

In summary, by incorporating nonlinear microion screening into a mean-field response theory of charged 
colloids in the primitive model, we have derived nonlinear corrections to the effective electrostatic inter- 
actions between hard spherical macroions in bulk colloidal suspensions. The key physical concept is that 
nonlinear screening entails both effective many-body interactions and essential corrections to the effective 
pair potential and the one-body volume energy. Effective triplet interactions are predicted to be always 
attractive, consistent with previous work |58l l8l) . The effective pair potential (r) , which in the linear 
(DLVO) theory has screened-Coulomb form, is shortened in range by the influence of nonlinear screening, 
but remains purely repulsive within the physically reasonable range of renormalized macroion charges. Pre- 

(2) I | 

dictions for tr ff (r) are in essentially perfect agreement with available ab initio simulation data |21| . The 
theory also predicts that triplet forces are negligible between weakly charged macroions, consistent with sim- 
ulation [2l]], but can be significant for higher macroion charges. Further simulations of more highly charged 
and concentrated macroions are now required to more severely test the theory. 

Analytical and numerical results confirm that nonlinear effects become qualitatively stronger with increas- 
ing macroion charge, increasing macroion concentration, and decreasing salt concentration. In the dilute 
limit of zero macroion concentration, but nonzero salt concentration, leading-order nonlinear corrections 
vanish. Perhaps the most practical application of the response theory, illustrated here, is to mapping out 
the parameter ranges within which linear theories can be trusted. Future work will explore implications of 
nonlinear screening for thermodynamic properties, in particular, the phenomenon of phase separation at low 
salt concentrations and the stability of deionized charged colloidal crystals. 
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APPENDIX: COMPARISON WITH RELATED THEORETICAL APPROACHES 



1. Response Theory vs. Poisson-Boltzmann Theory 

We demonstrate here that response theory, when combined with the random phase approximation (RPA) 
for the microion response functions (see Sec. IIII E|) is formally equivalent to Poisson-Boltzmann theory. The 
ensemble-averaged number density profile of positive microions, in the presence of the macroion potential 
0ext ( r ) ! is given exactly by the Euler-Lagrange equation [6^, H3 > 



P+( r ) 



Ai 



■ exp 



-#ze0ext(r) +cV" ) (r; [p+(r), p_(r)]) 



(A.l) 



which follows from minimization of the grand potential functional with respect to p+(r). In Eq. (|A.1|I . p + 
and A + denote the chemical potential and thermal de Broglie wavelength of the positive microions, (/> ex t(r) 

is the "external" electrostatic potential of the macroions, and <xt (r; [p+(r), p_(r)]) is the one-particle direct 
correlation function (DCF) of the positive microions, which is a functional of the inhomogeneous microion 

densities. Expanding c+\r; [p + (r), p_(r)]) in a functional Taylor series about the average (bulk) microion 
densities, n + and n_, we have 



"Ai 



■ exp 



Pze(f> cxt (r) + c ( | 1) (n + ,n_) + / dr' c^ 2 {(|r - r'|; n+, n_) [p+(r') - n+] 



dr' 



-,(2) 



|r - r'\; n + ,n-) [p-(r') - n_] + 



(A.2) 



(2) 

where (r; n + , n._), i,j = ±, are the bulk microion two-particle DCFs, which are related to the one-particle 
DCFs via 



q?(|r- r'|; n+,n_) 



lim 

p±( r )— »™± 



fe|"(r;kW.<'-(r)]) 
%(r') 



(A.3) 



We now make the mean-field random phase approximation |63| : (1) neglecting three-particle and higher- 
order correlations, i.e., truncating the series in Eq. (| A. 2|) . and (2) ignoring short-range pair correlations by 

f 21 

simply equating c\ ■ (r;rj + ,n_) to their asymptotic long-range limits, 



(2) 

Equation (|A. 21) then becomes 

p + (r) ~ n+exp -0 ^ze0 cxt (r) + J dr'u ++ (|r - r'|) [p+(r') - n+] 

+ 1 dr'«+_(|r-r'|)[p_(r')-n_]) 
= n + exp[— /3ze0(r)], 
where we have used n + = (e^ M + /A?) exp[c+ and have identified 

/?(■ 
dr' [p+(r') -n + - p_(r') + 

er-r 



(A.4) 



(A.5) 



(A.6) 



as the total electrostatic potential, due to both the macroions and the surrounding microions. Similarly, the 
density profile of negative microions is given by 



p_(r) 



n_ exp 



-0 ^-ze^tir) + J dr' U+ _(|r - r'|) [p+(r') - n+] 

+ J dr'u__(|r-r'|)[p_(r')-n_] 
n_ exp[fize<fi(r)]. 



(A.l) 
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Combining Eqs. (|A.5(1 and (|A.7(1 with the Poisson equation 



Air 

V </>(r) = {zep+{r) - zep_(r)} , 

e 



we obtain 



V 2 <b(r) 



Airze 



{n + exp[— /3ze0(r)] — n_ exp[/3ze</>(r)]} , 



(A.f 



(A.9) 



which is the PB equation for macroions in a symmetric z : z electrolyte. We conclude that the RPA-based 
response theory is formally equivalent to the mean-field PB theory. This is not surprising, given that both 
approaches neglect microion correlations. Response theory, however, provides a powerful framework for 
going beyond a mean-field description by systematically including microion correlations via more accurate 
approximations for the response functions of the microion plasma. 



2. Response Theory vs. Integral Equation Theory 



Here we show that linear response theory is equivalent to a linearized hypernetted-chain (HNC) approxima- 
tion in integral-equation theory. Substituting Eqs. I|72I74|I into the Fourier transform of Eq. J7SJ), the linear 
response expressions for the ensemble-averaged number density profiles of positive and negative microions 
are 



P+( r ) 
n+ 



Vm +(\ r - R *D + J dr ' i n + h ++(\ r - r 'D ~ n-h+-(\r - r'|)K l+ (|r' - R, 



(A.10) 



and 



Mr) 



= 1 + 0J2 v m+ (\r-R l \) + J dr' [n_/i — (|r — r'|) - n + h + -(\r - r'|)]u m +(|r' - R, ; 



On the other hand, from Eq. (|A.2|I . we have the exact relations 



p + (r) = n+exp 



-pJ2v m+ (\r-Ri\) + J dr'c ( f 2 l(|r-r'|)[ P+ (r')-n + ] 



+ / dr , c ( f 2 l(|r-r , |)[p_(r , )-n_] + --- 



and 



p_ (r) = n_ exp 



/3^> m+ (|r-R i |) + J dr' C L 2 i(|r-r'|)[p_(r')-n-] 



+ / dr'4 2 l(|r-r'|)[ P+ (r') -«+] + ■■■ 



(A.ll) 



(A.12) 



(A.13) 



Truncating the expansions on the right side of Eqs. (|A.12|I and l|A.13|l at the level of two-particle correlations 
amounts to the HNC approximation in integral-equation theory. If we further linearize the exponential 
functions (expanding and neglecting all but the first two terms) , substitute recursively for p + (r) and p_ (r) , 
and use the Ornstein-Zernike (OZ) relation for mixtures |63| 



MO^fM + E™* / dr' Cl (2) (|r-r'|)Mr'), 

7„ « 



(A.14) 



we recover Eqs. IjA.lOfl and i|A.ll|) . Thus, the linear response approximation is equivalent to a linearized- HNC 
closure of the OZ relation, while nonlinear response generates new closures. 



21 



[1] R. J. Hunter, Foundations of Colloid Science (Oxford University Press, Oxford, 1986). 

[2] P. N. Pusey, in Liquids, Freezing and Glass Transition, session 51, ed. J. -P. Hansen, D. Levesque, and J. Zinn- 

Justin (North-Holland, Amsterdam, 1991). 
[3] K. S. Schmitz, Macroions in Solution and Colloidal Suspension (VCH, New York, 1993). 
[4] F. Oosawa, Poly electrolytes (Dekker, New York, 1971). 
[5] Polyelectrolytes, ed. M. Hara (Dekker, New York, 1993). 

[6] B. V. R. Tata, M. Rajalakshmi, and A. K. Arora, Phys. Rev. Lett. 69, 3778 (1992). 
[7] K. Ito, H. Yoshida, and N. Ise, Science 263, 66 (1994). 

[8] H. Matsuoka, T. Harada, and H. Yamaoka, Langmuir 10, 4423(1994); H. Matsuoka, T. Harada, K. Kago, and 
H. Yamaoka, %bid 12, 5588 (1996); T. Harada, H. Matsuoka, T. Ikeda, and H. Yamaoka, ibid 15, 573 (1999). 

[9] F. Crohn and M. Antonietti, Macromolecules 33, 5938 (2000). 
[10] A. E. Larsen and D. G. Grier, Nature 385, 230 (1997). 
[11] G. M. Kepler and S. Fraden, Phys. Rev. Lett. 73, 356 (1994). 

[12] A. E. Larsen and D. G. Grier, Phys. Rev. Lett. 76, 3862 (1996); J. C. Crocker and D. G. Grier, Phys. Rev. Lett. 77, 
1897 (1996). 

[13] Y. Levin, Rep. Prog. Phys. 65, 1577 (2002). 

[14] C. N. Likos, Phys. Rep. 348, 267 (2001). 

[15] L. Belloni, J. Phys.: Condens. Matter 12, R549 (2000). 

[16] J.-P. Hansen and H. Lowen, Ann. Rev. Phys. Chem. 51, 209 (2000). 

[17] M. O. Robbins, K. Kremer, and G. G. Grest, J. Chem. Phys. 88, 3286 (1988). 

[18] E. J. Meijer and D. Frenkel, J. Chem. Phys. 94, 2269 (1991). 

[19] S. Auer and D. Frenkel, J. Phys.: Condens. Matter 14, 7667 (2002). 

[20] A.-P. Hynninen and M. Dijkstra, Phys. Rev. E 68, 21407 (2003). 

[21] R. Tehver, F. Ancilotto, F. Toigo, J. Koplik, and J. R. Banavar, Phys. Rev. E 59, R1335 (1999). 

[22] H. Lowen, P. A. Madden, and J.-P. Hansen, Phys. Rev. Lett. 68, 1081 (1992); H. Lowen, J.-P. Hansen, and 

P. A. Madden, J. Chem. Phys. 98, 3275 (1993); H. Lowen and G. Kramposthuber, Europhys. Lett. 23, 673 

(1993). 

[23] M. J. Stevens, M. L. Falk, and M. O. Robbins, J. Chem. Phys. 104, 5209 (1996). 

[24] P. Linse, J. Chem. Phys. 110, 3493 (1999); V. Lobaskin and P. Linse, J. Chem. Phys. Ill, 4300 (1999); P. Linse 
and V. Lobaskin, Phys. Rev. Lett. 83, 4208 (1999); P. Linse, J. Chem. Phys. 113, 4359 (2000); J. Rescic and 
P. Linse, J. Chem. Phys. 114, 10131 (2001); V. Lobaskin, A. Lyubartsev, and P. Linse, Phys. Rev. E 63, 20401 
(2001); V. Lobaskin and K. Qamhieh, J. Phys. Chem. B 107, 8022 (2003). 

[25] R. Messina, C. Holm, and K. Kremer, Phys. Rev. Lett. 85, 872 (2000); Eur. Phys. J. E 4, 363 (2001). 

[26] E. Allahyarov, I. DAmico, and H. Lowen, Phys. Rev. Lett. 81, 1334 (1998). 

[27] B. V. Derjaguin and L. Landau, Acta Physicochimica (USSR) 14, 633 (1941). 

[28] E. J. W. Verwey and J. T. G. Overbeek, Theory of the Stability of Lyophobic Colloids (Elsevier, Amsterdam, 
1948). 

[29] H. Graf and H. Lowen, Phys. Rev. E 57, 5744 (1998). 

[30] R. van Roij and J.-P. Hansen, Phys. Rev. Lett. 79, 3082 (1997). 

[31] R. van Roij, M. Dijkstra, and J.-P. Hansen, Phys. Rev. E 59, 2010 (1999). 

[32] M. J. Crimson and M. Silbert, Mol. Phys. 74, 397 (1991). 

[33] A. R. Denton, J. Phys.: Condens. Matter 11, 10061 (1999). 

[34] A. R. Denton, Phys. Rev. E 62, 3855 (2000). 

[35] P. B. Warren, J. Chem. Phys. 112, 4683 (2000). 

[36] D. Y. C. Chan, P. Linse, and S. N. Petris, Langmuir 17, 4202 (2001). 
[37] I. Rouzina and V. A. Bloomfield, J. Phys. Chem. 100, 9977 (1996). 

[38] B.-Y. Ha and A. J. Liu, Phys. Rev. Lett. 79, 1289 (1997); Phys. Rev. E 58, 6281 (1998). 
[39] M. J. Stevens Phys. Rev. Lett. 82, 101 (1999); Biophys. J. 80, 130 (2001). 

[40] B. I. Shklovskii, Phys. Rev. E 60, 5802 (1999); T. T. Nguyen, I. Rouzina, and B. I. Shklovskii, Phys. Rev. E 60, 
7032 (1999). 

[41] T. T. Nguyen, I. Rouzina, and B. I. Shklovskii, J. Chem. Phys. 112, 2562 (2000); T. T. Nguyen, A. Yu. Grosberg, 

and B. I. Shklovskii, J. Chem. Phys. 113, 1110 (2000). 
[42] W. M. Gelbart, R. F. Bruinsma, P. A. Pincus, and V. A. Parsegian, Physics Today 53 (Sept. 2000), p. 38. 
[43] G. N. Patey, J. Chem. Phys. 72, 5763 (1980). 
[44] L. Belloni, Phys. Rev. Lett. 57, 2026 (1986). 

[45] S. Khan and D. Ronis, Mol. Phys. 60, 637 (1987); S. Khan, T. L. Morton, and D. Ronis, Phys. Rev. A 35, 4295 
(1987). 

[46] M. D. Carbajal-Tinoco and P. Gonzalez-Mozuelos, J. Chem. Phys. 117, 2344 (2002). 



22 



L. B. Bhuiyan and C. W. Outhwaite, J. Chem. Phys. 116, 2650 (2002). 
S. N. Petris and D. Y. C. Chan, J. Chem. Phys. 116, 8588 (2002). 

J. A. Anta and S. Lago, J. Chem. Phys. 116, 10514 (2002); V. Morales, J. A. Anta, and S. Lago, Langmuir 19, 
475 (2003). 

H. H. von Griinberg, R. van Roij, and G. Klein Europhys. Lett. 55, 580 (2001). 

M. Deserno and H. H. von Griinberg, Phys. Rev. E 66, 011401 (2002). 

M. N. Tamashiro and H. Schiessel, J. Chem. Phys. 119, 1855 (2003). 

W. R. Bowen and A. O. Sharif, Nature 393, 663 (1998). 

J. J. Gray, B. Chiang, and R. T. Bonnecaze, Nature 402, 750 (1999). 

J. Dobnikar, R. Rzehak, and H. H. von Griinberg, Europhys. Lett. 61, 695 (2003); J. Dobnikar, Y. Chen, 
R. Rzehak, and H. H. von Griinberg, J. Phys.: Condens. Matter 15, S263 (2003). 

D. Goulding and J.-P. Hansen, Europhys. Lett. 46, 407 (1999). 

J.-P. Hansen, D. Goulding, and R. van Roij, J. Phys. IV France 10, Pr5-27 (2000). 

H. Lowen and E. Allahyarov, J. Phys.: Condens. Matter 10, 4147 (1998). 

M. E. Fisher, J. Stat. Phys. 75, 1 (1994); X. Li, Y. Levin, and M. E. Fisher, Europhys. Lett. 26, 683 (1994); 

M. E. Fisher, Y. Levin, and X. Li, J. Chem. Phys. 101, 2273 (1994). 

N. V. Sushkin and G. D. J. Phillies, J. Chem. Phys. 103, 4600 (1995). 

L. E. Gonzalez, D. J. Gonzalez, M. Silbert, and S. Baer, Mol. Phys. 99, 875 (2001). 

J. S. Rowlinson, Mol. Phys. 52, 567 (1984). 

J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, 2 nd ed. (Academic, London, 1986). 

N. W. Ashcroft and D. Stroud, Solid State Phys. 33, 1 (1978). 

J. Hafner, From Hamiltonians to Phase Diagrams (Springer, Berlin, 1987). 

N. W. Ashcroft, Phys. Lett. 23, 48 (1966). 

We note that expansion about zero macroion charge (it = 0) in Eq. I16|l is here the only natural choice, thus 
avoiding ambiguities encountered in other linearization schemes [Holl5lLl5^|. 

E. A. Allahyarov, L. I. Podloubny, and P. P. J. M. Schram, and S. A. Trigger, Physica A 220, 349 (1995). 

We note that the counterion-counterion pair correlation function h cc (r) in Eq. 14611 also appears in the distribution 
function approach of Chan et al. (cf. Eq. (45) of ref. |3fJ). 

Equation 1791 has also been derived by a density-functional approach in refs. |29l| and [3l|. 

I. S. Gradstein and I. M. Ryzhik, Table of Integrals, Series, and Products (Academic, New York, 1980). 

R. Klein, H. H. von Griinberg, C. Bechinger, M. Brunner, and V. Lobaskin, J. Phys.: Condens. Matter 14, 7631 
(2002). 

R. V. Durand and C. Franck, Phys. Rev. E 61, 6922 (2000). 

S. Alexander, P. M. Chaikin, P. Grant, G. J. Morales, and P. Pincus, J. Chem. Phys. 80, 5776 (1984). 

In ref. [Tof. a macroion valence of Z = 7300 is obtained by fitting the pair interaction between two isolated 

macroions with a DLVO potential. For bulk crystals, however, this valence somewhat exceeds the approximate 

charge-renormalization limit [74[. 

J. C. Neu, Phys. Rev. Lett. 82, 1072 (1999). 

J. E. Sader and D. Y. C. Chan, J. Coll. Int. Sci. 213, 268 (1999); Langmuir 16, 324 (2000). 

E. Trizac and J. L. Raimbault, Phys. Rev. E 60, 6530 (1999); E. Trizac, Phys. Rev. E 62, R1465 (2000). 

K. S. Schmitz, Phys. Chem. Chem. Phys. 1, 2109 (1999). 

R. P. Sear, Phys. Rev. E 62, 2501 (2000). 

C. Russ, H. H. von Griinberg, M. Dijkstra, R. van Roij, Phys. Rev. E 66, 011402 (2002). 

J. Z. Wu, D. Bratko, H. W. Blanch, and J. M. Prausnitz, J. Chem. Phys. 113, 3360 (2000). 

A. R. Denton (unpublished). 

C. J. Pethick, Phys. Rev. B 2, 1789 (1970). 

E. G. Brovman and G. Solt, Solid State Comm. 8, 903 (1970). 

S. P. Singh and W. H. Young, J. Phys. F: Metal Phys. 3, 1127 (1973). 

M. Rasolt and R. Taylor, Phys. Rev. B 11, 2717 (1975); L. Dagens, M. Rasolt, and R. Taylor, Phys. Rev. B 11, 
2726 (1975). 

A. A. Louis and N. W. Ashcroft, Phys. Rev. Lett. 81, 4456 (1998). 

From Eq. IA.1I on, to simplify notation, we dispense with angular brackets in denoting ensemble-averaged 
densities. 



23 




FIG. 2: Geometry for the physical interpretation of response theory (see Sec. lIIID|l . Vectors r and r' define center- 
to-center displacements of macroions. Vectors n, r2, and r-3 define points at which either the macroion "external" 
potential acts or a change in microion density is induced. 
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FIG. 3: Ensemble-averaged counterion density around a single macroion for macroion diameter a — 100 nm and 
valence (a) Z = 100, (b) Z = 500, at volume fraction ij — 0.01 and zero salt concentration. Dashed curves: linear 
response theory. Solid curves: nonlinear response theory (first-order correction). 
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FIG. 4: Effective pair interactions for macroion diameter a, valence Z, volume fraction 77, and salt concentration c s : 
(a) a = 100 nm, Z = 400, t] = 0.01, c s = 1 fj,M; (b) a = 100 nm, Z = 700, rj = 0.01, c s = 1 ^M; (c) a = 652 nm, 
Z = 5000, 77 = 0.0352, c s = 0.2 /xM (chosen to compare with ref. |Tcj|'). Dashed curves: linear response prediction. 
Solid curves: nonlinear response prediction. The insets show that for sufficiently weakly charged macroions (a) the 
effective pair interaction may be fit by a Yukawa potential (with a lower effective charge), while for more highly 
charged macroions (b,c) deviations from linear behavior can be significant. 
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FIG. 5: Total interaction potential energy for two macroions, of diameter a — 106 nm and valence (a) Z = 200, 
(b) Z = 700, in a cubic box of length 530 nm with periodic boundary conditions at zero salt concentration. The 
potentials are shifted to zero at maximum macroion separation. Dashed curves: linear response prediction. Solid 
curves: nonlinear response prediction. Symbols: ab initio simulation data |2'lll . 
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FIG. 6: Map of nonlinear deviations from linear response theory for macroions of diameter a = 100 nm and 
valences, from top to bottom, Z — 500, 600, 700, for fee crystal structure. Systems with macroion volume fractions 
rj and salt concentrations c s above the respective curves deviate from the linear response pair potential at the fee 
nearest-neighbor distance by at least (a) 1 fcsT or (b) 0.1 ksT. 
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FIG. 7: Effective three-body interaction between three macroions, arranged in an equilateral triangle of side length 
r, with macroion diameter a = 100 nm, valence Z — 500 (open circles), Z = 700 (filled circles), volume fraction 
r\ = 0.01, and salt concentration c s — 1 fiM. Computed by Monte Carlo integration of Eq. I|82|l . with numerical errors 
comparable to symbol size. 
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FIG. 8: Predicted force on a macroion in an equilateral-triangle arrangement of three macroions, each of diameter 
a = 106 nm and valence (a) Z = 700, (b) Z = 1000, in a cubic box of length 1000 nm with periodic boundary 
conditions at zero salt concentration. Dashed curves: sum of linear response effective pair forces. Solid curves: sum 
of nonlinear effective pair forces. The symbols are theoretical predictions for the effective triplet force (open circles) 
and the total (pair plus triplet) effective force (filled circles). 



